COMMENT ⊗   VALID 00007 PAGES
C REC  PAGE   DESCRIPTION
C00001 00001
C00002 00002	ENTRY incr_prime
C00003 00003	∂ Declarations for SIEVE.
C00004 00004	INTERNAL INTEGER PROCEDURE incr_prime(
C00007 00005	    ∂ Go looking for primes.
C00011 00006	    ∂ Got all the primes needed, or else ran out of table.
C00012 00007	END "PRIMER".
C00013 ENDMK
C⊗;
ENTRY incr_prime;
BEGIN "PRIMER"
REQUIRE "HEADER.SAI[REV,KS]" SOURCE_FILE;

∂ Ken Shoemake.  June 1977.
This module finds the k'th prime above or below a given integer, or the nearest
prime if k is zero.
;
∂ Declarations for SIEVE.;
∂ These things currently found in PSIEVE.SAI[REV,KS];
EXTERNAL INTEGER #SIEVE;
EXTERNAL INTEGER ARRAY SIEVE[0:(#SIEVE DIV 36)]; ∂ Bounds are ignored.;

INTERNAL INTEGER PROCEDURE incr_prime(
	INTEGER n, incr(0));
∂ Increments or decrements n by incr prime numbers.  If incr is 0, returns
nearest prime to n.  If n is outside the bounds of the prime table, it is
returned as the result, and an error message is printed.  Otherwise the
result is guaranteed to be prime.  An error message will be printed if the
attempt to increment/decrement n exceeds the range of the table, and the
maximum/minimum prime in the table will be returned as the result.  The
table is stored as a bit vector, 36 bits to a word, where a one bit means
that index is prime.  Both the bit vector (SIEVE) and the maximum prime
in it (#SIEVE) are defined EXTERNALly in PSIEVE.REL, which is produced from
PSIEVE.SAI, which is generated by executing PSGEN.SAI.  Counting the high
order bit in a word as bit 0 and the low order bit as bit 35, the bit
corresponding to a given integer n will be bit (n MOD 36) in SIEVE[(n DIV 36)].
;
    BEGIN "incr prime"
    INTEGER wrd, bit, dir;

    IF NOT (2 ≤ n ≤ #SIEVE)
      THEN BEGIN
	PRINT(↓,"incr_prime: 2 ≤ n ≤ ",#SIEVE,↓);
	RETURN(n);
	END;

    wrd ← n DIV 36;
    bit ← n MOD 36;

    IF incr = 0
      THEN BEGIN "nearest prime"
	INTEGER top, bot;
	IF (SIEVE[wrd] LSH bit) < 0	∂ Make n's bit the sign bit and test it;
	  THEN RETURN(n);
	top ← incr_prime(n,1);
	bot ← incr_prime(n,-1);
	IF ((top+bot) DIV 2) ≥ n	∂ Choose the nearer prime, lower if tie;
	  THEN RETURN(bot)
	  ELSE RETURN(top);
	END "nearest prime";

    ∂ Go looking for primes.;
    IF incr < 0
      THEN BEGIN
	dir ← -1;
	incr ← -incr;
	END
      ELSE BEGIN
	dir ← 0;
	wrd ← wrd LOR ((wrd-((#SIEVE DIV 36)+1)) LSH 18);
	END;


	START_CODE
	DEFINE inc=0, tp=1, ac=2, bt=ac+1, wd=4, dr=5, sv=6;
	LABEL Xcrement,Decr,Mask,Nope,GotOne,Pinpoint,GotIt,IncrWd,Lose,TheEnd;

		MOVEI	sv,ACCESS(SIEVE[0]);	∂ Need to have this address;
		TLO	sv,wd;		∂ Will be indexed by wd;
		MOVE	inc,incr;	∂ These should come →after← ACCESS ;
		MOVE	bt,bit;		∂   to avoid getting clobbered;
		MOVE	wd,wrd;
		MOVE	dr,dir;

	    Xcrement:MOVN    tp,bt;	∂ Go up/down looking for another prime;
		    JUMPL   dr,Decr;	∂ Mask off bits behind current position;
		    SETO    ac,;
		    LSH     ac,-1(tp);	∂ Sets ac←00...0x1...11, x=current,set to 0;
		    JRST    Mask;
	    Decr:   MOVSI   ac,'400000;
		    ASH     ac,1(tp);	∂ Sets ac←11...1x0...00, x=current,set to 0;
	    Mask:   AND     ac,@sv;	∂ Address is effectively SIEVE[wd];
		    JUMPN   ac,GotOne;	∂ Is there a prime in the same word?;

		Nope:   XCT     IncrWd(dr);	∂ No primes yet, try next word;
			SKIPN   ac,@sv;	∂ Any primes in it (SIEVE[wd]) ?;
			JRST    Nope;

	    GotOne: JUMPGE  dr,Pinpoint;∂ Got a prime, but exactly what bit is it?;
		    MOVN    tp,ac;	∂ Get the proper one for this direction;
		    ANDM    tp,ac;	∂ Pick out lowest one bit, else highest;
	    Pinpoint:JFFO   ac,GotIt;	∂ Put the bit number in bt (=ac+1);
	    GotIt:  SOJG    inc,Xcrement;	∂ Go back for more, unless sated;

		HRRZM	wd,wrd;		∂ Return word and bit where found prime;
		MOVEM	bt,bit;
		MOVEM	inc,incr;	∂ Will be zero, indicating total success;
		JRST	TheEnd;

	∂ The following instructions are not executed in-line, but by XCT at Nope;
		SOJL    wd,Lose;	∂ Decrement and test against 0;
	IncrWd: AOBJP   wd,Lose;	∂ Increment and test against #SIEVE/36;

	Lose:	HRRZM	wd,wrd;		∂ Return word where lost;
		MOVEM	inc,incr;	∂ Will be non-zero, indicating failure;

	TheEnd:
	END;
    ∂ Got all the primes needed, or else ran out of table.;
    ∂ Number of primes left to find will be in incr.;

    IF incr = 0
      THEN RETURN(bit+36*wrd)
      ELSE BEGIN
	PRINT(↓,"incr_prime: 2 ≤ n ≤ ",#SIEVE,↓);
	IF dir = -1		∂ If going down, then minimum prime, else maximum;
	  THEN RETURN(2)
	  ELSE RETURN(#SIEVE);	∂ #SIEVE is the largest prime in table. Sorry;
	END;

END "incr prime";
END "PRIMER".